function D=correct_planar_grads(filename, fixed_tra_fname,factor)

% correct_planar_grads(filename)
%
% - corrects grad.tra so that Neuromag-306 planar gradiometers are now
% stored with the correct magnitiude
%
% MWW 


%% try loading in file

try
    D = spm_eeg_load(filename);
catch
    warning('Unable to run correct_planar_grads on file %s; probably not an MEEG file',filename);
    return;
end

if sum(ismember(unique(D.chantype),'MEGGRAD')) % HL Mod 1.3 - Prevents CTF data being corrected.
    error('Axial Gradiometers Detected! Are you using CTF data? If so you should not be correction gradiometers.')
end

%% correct gradiometer definition;

load(fixed_tra_fname);

%1. correct D.sensors
grad = D.sensors('MEG');

% check to see if data has already been montaged - if it has then error
if sum(abs(grad.tra(1,:))>0) > 2,
    error('correct_planar_grads: Can not do grad baseline correction on data that has already been montaged.');
end;

x=factor*fixed_tra;
tmp1=[3:3:size(x,1)];
tmp2=[5:5:size(x,2)];

for i=1:length(tmp1),x(tmp1(i),tmp2(i))=1; end;

grad.tra=x;

%%%%
x=factor*fixed_tra;
tmp1=[3:3:size(x,1)];
tmp2=[5:5:size(x,2)];

for i=1:length(tmp1),x(tmp1(i),tmp2(i))=1; end;

grad_opp=grad;
grad_opp.tra=x;
grad_opp.tra(1,1:2);

grad.tra_baseline_corrected=1;

%%%%

D = sensors(D,'MEG',grad);

%2. correct forward models
if isfield(D,'inv');
    for i = 1:length(D.inv),
        %2.1 correct sensors in datareg
        if isfield(D.inv{i},'datareg')
            for j = 1:length(D.inv{i}.datareg),
                if strcmp(D.inv{i}.datareg(j).modality,'MEG')
                    disp(['Before correction: tra=']);
                    disp(D.inv{i}.datareg(j).sensors.tra(1:5,1:5));

                    D.inv{i}.datareg(j).sensors.tra = grad.tra;                                        
                    D.inv{i}.datareg(j).grad_opp=grad_opp;
                    
                    disp(['After correction: tra=']);
                    disp(D.inv{i}.datareg(j).sensors.tra(1:5,1:5));
                end
            end
        end
        
        %2.2 remove previously computed leadfield so that it is recomputed
        %    in future inversions
        if isfield(D.inv{i},'gainmat')
            fname = fullfile(D.path, D.inv{i}.gainmat);
            if exist(fname,'file') %delete the lead field file
                cmd = ['rm ' fname];
                unix(cmd);
            end
            D.inv{i} = rmfield(D.inv{i},'gainmat'); %remove the reference to it
        end
    end
end

save(D);
fprintf('Successfully corrected %s \n',filename);
                    


    